
#include "glodef.h"
#include "errlog.h"
#include "fe.h"
#include "fe_matr.h"

#include <stdlib.h>

#define LAPLACE_ID           0x4CBF064C

typedef struct
{
    unsigned id;
    unsigned eltype;
    grid3d_t const *grid;
    fe_matr_p laplace;
} laplace_t;

static laplace_t *get_ldata(void *data)
{
    laplace_t *l=data;
    if((data==NULL)||(l->id!=LAPLACE_ID))
        error_msg("laplace.c: invalid laplace_t structure passed.\n");
    return l;
}

void SolveLaplace_scalar(void *ldata, scalar U, scalar F, scalar T,
                         unsigned F_tested, unsigned reassemble)
{
    laplace_t *ld=get_ldata(ldata);
    if(!F_tested)
    {
        test_scalar(ld->grid, ld->eltype, F, T);
        copy_scal(T, F);
    }
    if(reassemble)
        fe_matr_assemble(ld->grid, 1, &ld->laplace);
    PrintLog("Solving Laplace problem.\n");
    fe_matr_solve_scalar(ld->laplace, U, F, T, 0.0);
    return;
}

void SolveLaplace_vector(void *ldata, vector U, vector F, vector T,
                         unsigned F_tested, unsigned reassemble)
{
    laplace_t *ld=get_ldata(ldata);
    if(!F_tested)
    {
        test_vector(ld->grid, ld->eltype, F, T);
        copy_vec(T, F);
    }
    if(reassemble)
        fe_matr_assemble(ld->grid, 1, &ld->laplace);
    PrintLog("Solving Laplace problem.\n");
    fe_matr_solve_vector(ld->laplace, U, F, T, 0.0);
    return;
}

void *InitLaplace(grid3d_t *grid, unsigned eltype)
{
    LOG_MPI_WTIME_START;

    laplace_t *ld;
    InitFE(1, &eltype, grid);

    PrintLog("Initialize Laplace solver (%s) begin.\n", __func__);

    ld=malloc(sizeof(laplace_t));
    if(ld==NULL)
        error_msg("laplace.c: memory allocation error.\n");
    ld->id=LAPLACE_ID;
    ld->eltype=eltype;
    ld->grid=grid;
    ld->laplace=fe_matr_alloc(eltype, local_stif, eo_bc_dir, grid, TRUE,
#if defined(USE_CSC)
            "ss"
#else
            NULL
#endif
            );
    // fe_matr_structure(ld->grid, 1, &ld->laplace);
    fe_matr_assemble(ld->grid, 1, &ld->laplace);

    LOG_MPI_WTIME_STOP;

    return ld;
}

void ReassembleLaplace(void *ldata)
{
    laplace_t *ld=get_ldata(ldata);
    fe_matr_assemble(ld->grid, 1, &ld->laplace);
}

void DoneLaplace(void *ldata)
{
    laplace_t *ld=get_ldata(ldata);
    fe_matr_free(ld->laplace);
    free(ld);
    PrintLog("Laplace solver done.\n");
}


